•O-A042  728  WISCONSIN  UNIV  MADISON  MATHEMATICS  RESEARCH  CENTER  F/6  3/9 

HuN^77^”S  CERTAIN  THREE-OI-ETC (U) 

UNCLASSIFIED  ® DAA6S9-75-C-002. 


DOC  FILE  COPY 


MRC  Technical  Summary  Report  # 1752 

USING  THE  METHOD  OF  ORTHOGONAL 
COLLOCATION  FOR  CERTAIN  THREE- 
DIMENSIONAL  PROBLEMS  OF  STELLAR 
STRUCTURE  ? 

M.  J.  Miketinac  and  S.  V.  Parter 


Mathematics  Research  Center  / ' 
University  of  Wisconsin-Madisg’n  / 
610  Walnut  Street  > j 

Madison,  Wisconsin  5370^6  7 


1 /// 
/ I 


June  1977 


Received  January  14,  1977 


i >■ 


O P, 


Sponsored  by 

U.S.  Army  Research  Office 
P.O.  Box  12211 
Research  Triangle  Park 
North  Carolina  27709 


Approved  for  public  release 
Distribution  unlimited 


Office  of  Naval  Research 
Arlington,  Virginia  22217 


f 


UNIVERSITY  OF  WISCONSIN  - MADISON 
MATHEMATICS  RESEARCH  CENTER 


USING  THE  METHOD  OF  ORTHOGONAL  COLLOCATION 
FOR  CERTAIN  THREE-DIMENSIONAL  PROBLEMS  OF  STELLAR  STRUCTURE? 

t 

M.  J.  Miketinac  and  S.  V.  Parter 


Technical  Summary  Report  #1752 
June  1977 


ABSTRACT  | 

The  method  is  developed  for  two  specific  problems:  (i)  computation  of 

I 

the  structure  of  the  primary  component  (assumed  to  consist  of  a polytropic  gas) 

in  a synchronous  close  binary  system  and  (ii)  search  for  non-axisymmetric  con-  i' 

figurations  of  differentially  rotating  polytropes.  In  both  cases  the  structure  v 

i' 

equations  reduce  to  a mildly  non-linear  elliptic  partial  differential  equation  in  i 

three  dimensions  with  boundary  conditions  at  the  center,  on  a sphere  containing 
the  star  and  involving  a 'free'  boundary.  The  present  method  has  several  advan-  j 

tages  over  the  'standard'  methods  (namely,  improvements  of  Chandrasekhar's  ^ 

perturbation  analysis).  The  most  important  of  these  are  consistency  and  easier  ^ 

application  to  real  stars.  However,  the  method  becomes  computationally  inef- 
ficient when  used  for  computing  of  configurations  with  strong  angular  dependence. 

In  such  cases  (related)  Galerkin  methods  offer  significant  advantages. 
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USING  THE  METHOD  OF  ORTHOGONAL  COLLOCATION 
FOR  CERTAIN  THREE-DIMENSIONAL  PROBLEMS  OF  STELLAR  STRUCTURE  ? 

M.  J.  Miketinac^and  S.  V.  Barter 

1.  Introduction 

Several  successful  attempts  at  numerical  solution  of  stellar  structure  problems  have  been  made 
in  the  past  dozen  years  (see  Strittmatter,  1969  and  Papaloizou  and  Whelan,  197  3).  Perhaps  the  most 
elegant  of  these  methods  has  been  the  one  devised  by  Stoeckly  (1965).  Most  of  these  methods  suffer 
from  serious  limitations  on  the  range  of  their  applicability,  or  undue  complexity  and/or  unproven  con- 
vergence properties.  Even  so,  none  of  them  has  ever  been  applied  to  three-dimensional  problems. 

In  this  Investigation,  Stoeckly's  formulation  of  the  structure  problem  for  rotating  stars  is 
altered  to  include  most  of  the  other  stellar  structure  problems  (section  2).  Then,  his  numerical  method, 
which  is  in  fact  a special  form  of  the  method  of  orthogonal  collocation  (used  extensively  in  theoretical 
chemistry:  Finlayson,  1972),  is  improved  and  generalized  to  three  dimensions  (sections  3 and  .1).  The 
method  is  developed  for  polytropic  models  of  stars  and  in  its  present  form  it  is  applicable  only  to  poly- 
tropes  with  the  polytropic  index  n>  1.  However,  it  would  be  fairly  straightforward  to  modify  the  meth- 
od in  such  a way  that  polytropes  with  n < 1 and,  also,  matter  more  complicated  than  polytropes 
(Miketinac,  1976)  could  be  treated.  Convergence  of  the  method  is  virtually  assured  through  the  work  of 
Vainikko  on  perturbed  Galerkin  methods  (1972,  1967).  The  method  produces  results  (section  5)  in  satis- 
factory agreement  with  the  known  results  for  uniformly  rotating  configurations  of  polytropes  (which  is 
actually  a two-dimensional  problem).  However,  the  method  has  an  undesirable  feature  making  its  use 
on  the  computer  much  more  expensive  than  the  'perturbed'  Galerkin  method.  This  last  method  is,  there- 
fore, the  recommended  method  for  computations  of  three-dimensional  stellar  structure  models.  It  is 
shown  in  section  3 that  the  method  of  orthogonal  collocation  is  in  fact  an  equivalent  - but  computation- 
ally less  efficient  (shown  in  section  4)  - form  of  the  'perturbed'  Galerkin  method.  It  is  likely  tliat  the 
'perturbed'  Galerkin  method  as  formulated  in  sections  3 and  4 will  have  applications,  also,  outside 
the  stellar  structure  problems. 

This  report  has  also  been  issued  by  the  Computer  Sciences  Department,  University  of  Wiscon- 
sin as  T. R.  No.  298. 

t On  leave  from  Department  of  Applied  Mathematics.  University  of  Cape  Town. 

“Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-7 5-C-0024  and  by  the  O.  N.  R.  under 
Contract  No.  N00014-76-C-0341. 
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2,  Structure  Equations  for  a Polytropic  Model 

It  can  be  shown  that  equations  describing  the  configuration  of  a self-gravitating,  polytropic 
fluid  modelling  a star  in  equilibrium  under  the  influence  of  some  disturbing  force  reduce  to 

= -e"  , (2.1) 

e = V - + 1 + . (2.2) 

All  symbols  in  these  equations  - called  the  structure  equations  - represent  dimensionless  quantities 
The  equations  are  formulated  in  the  spherical  polar  coordinate  system  with  the  center  coin- 

ciding with  the  center  of  symmetry  (usually  the  center  of  mass)  and  with  the  axes  oriented  along  the 
lines  of  symmetry  of  the  star.  In  this  coordinate  system  the  Laplace  operator,  , is 


j.  J_r  -1/1  2,  0 ^ 

3^2  X 0x  ^2f  ^ 


where  p = cos^  . Symmetries  of  the  star  model  depend  on  the  given  function 


I 


= Vp(\,x,p,v>)  (2.4) 

which  is  the  potential  describing  the  disturbing  force.  It  is  assumed  that  the  function  vanishes 
at  X = 0,  and  that  the  positive,  arbitrary  parameter  X.  - characterizing  the  strength  of  disturbance- 
is  such  that  for  X = 0,  Vj^=  0 (in  this  case  the  star  is  undisturbed  and  assumes  a spherical  shape). 
The  function  -V(x,p,</>)  is  essentially  the  gravitational  potential  of  the  model  and  = V(x=0,  p,  e). 
The  parameter  n - called  the  polytropic  index  - is  arbitrary  in  the  range  0<  n < s (it  specifies  the 
chemical  composition  of  the  model)  The  equation  (2.2)  - known  as  the  equation  of  hydrostatic  sup- 
port - Is  valid  only  inside  (defined  by  0 > 0)  the  star;  the  star's  surface  Is  defined  as  that  surface  on 
which  0(x,p,  ^)  = 0 . Since  0 is  not  known  in  advance,  this  gives  the  structure  problem  a 'free 
boundary'  aspect.  The  function  0 , where  positive,  is  related  to  the  density  and  is  such  that 

0(x  = 0,p,<p)  = 1 . (2.  5) 

For  convenience,  (2.  2)  will  be  taken  to  define  0 throughout,  but  then  in  equation  (2.1)  0 can  only 
be  positive,  i.e.  (2.1)  is  replaced  by 

v\  = - e"  , (2.1') 

where 
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when  0 > 0 
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(2.6) 


0 when  0 < 0 

From  this  it  follows  that  (k+2)-th  derivative  of  the  potential  V with  respect  to  or  v>  exists  at 

the  star's  surface  only  if  k < n,  otherwise  (k+2)-th  derivative  is  singular  there. 


The  structure  equations  must  be  solved  under  the  following  conditions:  (i)  at  the  center. 

X - 0,  VV  is  continuous,  (ii)  at  the  surface  of  the  star  VV  is  continuous  and  (iii)  at  the  outside 
VV-  0 as  x-w.  A solution  V^,V,0  of  the  structure  equations  (2.1')and  (2.  2)  and  the  above 
boundary  conditions  is  a possible  configuration  of  the  model  and  is  characterized  by  a specific  value 
of  (n,\).  All  quantities  of  astrophysical  interest  can  be  computed  from  this  solution.  The  most  in- 
teresting values  of  n are  n = 1.  5 (neutron  stars,  Cowling  model)  and  n = 3.  0 ( supermassive  stars). 

Examples  of  astrophysical  systems  whose  structure  can  be  modelled  with  equations  (2.)'). 

(2.  2)  and  the  associated  boundary  conditions  include  (i)  (uniformly  or  differentially)  rotating  stars, 

(ii)  magnetic  stars  (with  or  without  rotation),  and  (iii)  binary  stars.  In  the  case  of  uniform  rotation, 
the  function  is  given  by 

Vp  = 2 x^  (1-p^)  (2.7) 

where  the  parameter  \ is  now  the  (dimensionless)  angular  velocity  . Differential  (oi  non-uniform) 
rotation  can  be  described  by  (Stookly,  1965) 


u - f I x^l  - ^i‘’-) 

''d  - lb  f ® J ' <2-8) 

where  the  constant  b is  a parameter  of  non -uniformity  of  rotation,  ranging  from  b = 0,  for  uniform 
rotation,  to  b — 1,  which  approximates  spatial  dependence  of  the  centrifugal  potential  possibly 
arising  during  contraction  from  a uniformly  rotating  mass  of  initially  homogeneous  density.  Auchmuty 
and  Beals  (1971a,  1971b)  consider  the  question  of  existence  and  regularity  of  solutions  for  some  models 
of  rotating  stars  and  show,  using  the  variational  method,  that  physically  reasonable  solutions  exist 
In  most  cases.  However,  their  conditions  arc  so  stringent  that  the  (physically  interesting)  polytrope 
with  n =3.0  is  excluded  and  the  uniform  rotation  is  inadmissible.  Auchmuty  (1974,  1975)  formulates 
a numerical  scheme  based  on  the  variational  method  for  obtaining  solutions  of  the  structure  problem 
for  rotating  stars.  He  shows  that  the  method  converges,  but  results  of  his  computations  and  compari- 
sons with  other  methods  have  yet  to  appear. 
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For  the  primary  in  a synchronous  binary  system  there  are  two  versions  of  the  disturbing  poten- 


tial. In  the  so  called  first-order  theory  (Chandrasekhar,  193J) 

4 


, mq  V , X J r,  , . . > . m(l+q)  , , X.2  2 . 

'd  = 4~^X  ,^/x'  7sin4cos».)+  -)  sin  «, 

J — fa 


(2.9) 


where  X is  the  separation  between  the  centers  of  mass,  m is  the  mass  of  the  undisturbed  primary 

(i.e.  the  mass  when  ^ = 0)  and  q = — , m being  the  mass  of  the  secondary.  In  the  second-order 
X m 


theory  (Martin,  1970) 


,,  mq  V , xj  „ , , . , . m(l+q)  1 , X ,2  . 1 

''d  = ^ < X>  7®*""  ^ ^ ‘ X ' 

J — fa 


(2.10) 


3 q(b+c-2a)  x „ , 

- - — — - Pj(sin^  cosv>), 

4irX 

where  the  new  symbols  a,c,b  are  the  moments  of  inertia  of  the  primary  about  the  axis  pointing  to  tlic 
secondary,  the  rotation  axis  and  the  third  axis,  respectively.  All  physical  quantities  in  (2.9)  and 
(2. 10)  are  dimensionless,  and,  m and  (b+c-2a)  in  (2.10)  are  obtained  from  the  first-order  theory. 
Expressions  (2.9)  and  (2.10)  are  approximate  models  of  the  real  system  consisting  of  two  deformable, 
extended  bodies  interacting  gravitationally  while  each  body  rotates  uniformly  at  the  same  rate  as  it 
revolves:  the  axes  of  rotation  and  revolution  being  all  parallel.  The  first  order  expression  (2.9)  in- 
cludes all  effects  of  the  secondary  on  structure  of  the  primary  up  to  the  order  of  magnitude  0(  v^), 

Xq 

where  v = with  x^  being  the  radius  of  the  undistorted  primary.  The  second-order  expression 

(2.10)  includes  effects  of  magnitude  up  to  O(v^).  To  these  orders  of  magnitude  structure  of  the 
primary  is  independent  of  the  details  of  structure  of  the  secondary. 

Both,  Chandrasekhar  and  Martin,  combine  (2.1)  and  (2.  2)  into  an  equation  for  6,  and  then 
seek  a solution  in  the  form  of  a perturbation  expansion 

J ,, 

0(x,p,«>)  = 0 (X)  + ^ V e (x,yi,<p) 
k=3 

where  the  cutoff  J is  equal  to  9 in  the  first-order  theory,  7 in  the  second-order  theory.  This  ex- 
pansion can  not  bo  uniformly  valid  throughout  the  primary,  because  O^fx)  = 0 at  x = x^  (the  so 
called  Emden  radius)  inside  the  primary.  Martin  (1970)  argues,  however,  that  the  expansion  with 

g 

J = 7 is  valid  around  x = x^  If  the  polytropic  index  n > j , so  that  his  results  for  the  polytrope 
n = 3.0  should  be  quite  accurate.  But  the  same  argument  shows  that  even  the  J = expansion  is 
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valid  around  x - only  if  n > --  = 2 leaving  out  tlie  other  interesting  polytiope  (for  which  n-l. 
This  inability  of  the  (regular)  peiturbation  analysis  (directly  traceable  to  the  non-unifonnii^'  oi 
above  expansion)  to  cope  with  both  n = 1.  5 and  n = 3.0  polytropes  suggests  that  solution  o;  '.hi- 
structure  problem  should  be  sought  either  numerically  or  using  techniques  of  th.e  singular  perturnatu' 
analysis.  The  first  alternative  provides  the  major  motivation  for  examining  the  method  of  orthoao  ud 
collocation  in  this  paper.  For  the  second  alternative  in  the  case  of  uniform  rotation  see  Smiti'. 

1976). 

In  all  cases,  (2.7)  to  (2.10)  the  following  operations 


tp  2rt-v 


(2.  11) 


do  not  change  the  system  as  a consequence  of  symmetries  built  in  the  expressions  for  the  dist  .rci:;c 
potential.  This  means  that  a (numerical)  solution  of  the  structure  problem  need  to  bo  found  on!v  r- 
one  half  of  the  space  above  the  equatorial  plane.  A further  restriction  of  the  domain,  in  which  i 
solution  is  sought,  can  be  obtained  by  a reformulation  of  the  outside  boundary'  conditio:'. 


For  every  physically  interesting  solution  (such  solutions  are  assumed  to  exist)  of  the  struc- 
ture problem,  there  is  a sphere  of  some  finite  radius  completely  containing  the  region  where  the  func- 
tion G is  positive.  Let  x^  be  the  greatest  of  these  radii  for  a certain  sequence  of  configurations . 
Then  outside  this  sphere  the  potential  satisfies  the  Laplace's  equation  and  it  is  expansiule  in  the 
form 


where 


V(X,p,vr)  = Yj 


t 

y 

f =0  m=0 


V*  (X) 

f m 


^ ( p.  t tp ) 

f m 


(2.12) 


when  m ^ 0 and  for  m = 0 


y;;*(p,.p) 


j 2f +1 

( f-m)'. 

y 2ir 

(f  +m)'. 

/ 2f+l 

/ 4n 

P j (p)  cos  m (p 


and  the  summation  in  (2.12)  is  such  that  m+f  is  always  even.  This  restriction  and  the  fact  that  only 
cosines  of  the  angle  <p  appear  are  due  to  symmetries  (2. 11).  The  coefacients  satisfy  the  following 
differential  equation 


r 2d  <(<+!)  , X' 

There  is  only  one  solution  which  remains  finite  as  x-*">o;  this  solution  is 

, > -(f+1) 

V.  (x)  a X 
f m 

and  on  the  sphere  of  radius,  x , completely  outside  the  star  it  can  be  expressed  as 

5;^  *-;:r ''-m'V  ■ "• 

H 


(2.13) 


For  fixed  n the  solution  of  system  of  equations  (2.1'),  (2.2),  and  (2.13)  depends  only  on  the 
parameter  \ . This  sequence  of  configurations  'starts'  with  X = 0,  the  spherically  symmetnc  con- 
figuration, and  'ends'  with  the  critical  configuration,  when  X = X^.  For  n >1  (which  includes  the 
two  interesting  cases)  the  critical  configuration  for  uniformly  rotating  stars  (James,  1964)  and  syn- 
chronous binary  systems  (Martin,  1970)  is  characterized  by  the  fact  that  the  sum  of  all  forces  at  the 

equator  just  balances  to  zero,  i.e. 

80 


-(x^,0,0)  = 0. 


(2.  14) 

Locating  the  critical  configuration  is  the  main  purpose  of  this  investigation.  It  will  bo  done  by 
chec'niny  the  left-hand  side  of  (2.14)  for  every  computed  configuration  of  the  sequence  characterized 
by  fixed  n and  adjusting  the  value  of  X in  such  a way  as  to  make  it  as  close  to  zero  as  possible 
but  keeping  it  negative. 

In  the  case  of  non-uniform  rotation,  when  the  disturbing  potential  is  given  by  (2.8),  the  se- 
quence of  axisymmetric  equilibrium  configurations  for  polytropes  with  n >1  does  not  necessarily  end 
with  a configuration  in  which  the  effective  gravity  is  zero  at  the  equator  (i.e.  when  the  condition 
(2.11)  is  satisfied).  Stoeckly  (1965)  provides  some  evidence  that  for  n = 1.  5 and  b close  to  1 a 
point  of  bifurcation  is  reached  and  the  sequence  continues  with  non-axisymmetric  configurations. 
James  (1964)  has  shown  numerically  that  such  bifurcation  points  exist  for  uniformly  rotating  polytropes 
with  the  fiolytropic  index  n < 0.808.  The  Maclaurin  and  Jacobi  ellipsoids  (which  can  be  thought  of 
as  uniformlv  lotaMng  polytropcs  with  n = 0)  are  a famous  example  of  bifurcation.  A numerical  search 
tor  f iturcaticin  iioinis  and  continuatioji  of  the  sequence  of  solutions  along  the  nev,'  branches  could  be 
mounted  usine  l.('  ideas  of  H.  13.  Keller  (1976).  This  provides  the  second  major  motivation  for  looking 
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at  the  method  of  orthogonal  collocation  as  a means  of  computing  three-dimensional  stellar  configura- 
tions. 

^ 3.  The  Method  of  Orthogonal  Collocation  as  a 'Perturbed'  Galerkin  Method 

I The  sequence  of  solutions  of  the  system  of  non-linear  equations  (2.1'),  ( 2.  2)  and  ( 2. 1 .3) . 

i 

I characterized  by  fixed  value  of  n and  variable  value  of  X.  , can  be  generated  from  the  spherically 

: symmetric  solution,  X = 0,  (which  can  always  be  computed  as  an  initial-value  problem  for  a sinclc 

! , 

ordinary  differential  equation  of  second  order,  Miketinac,  1976)  by  using  a combination  of  the 
Newton's  linearization  and  the  trial  free-boundary  method  (Cryer,  1976).  This  is  possible  because  a 
■*  guess  for  the  solution  V^,V  provides,  because  of  (2,2),  also,  a 'guess'  for  the  free  boundary’.  The 

procedure  consists  of  two  steps 

■I  (i)  approximating  an  unknown  solution  X,  V ,V  by  X,  V ,V  where  V ,'2  is  a know.i  solution  with 

c c c 

^ X = X - 6X  and  6X  is  'small', 

(ii)  improving  this  approximation  by  X,  V + 6V  , V+  6V  whore 
. c c 

, V^6V  + 60"  = -V^V  - o"  . 

' + + 

The  second  stop  is  iterated  with  V ,V  replaced  by  V + 6V  , V + f.V  until  certain  convc  r.conces 

c c c 

criteria  are  satisfied.  In  the  second  step 

G = V + 1 - + Vj^(X)  (3.1) 

(which  defines  and  'moves'  the  boundary),  and 

60"  = n 6”“’  (6V  - 6V  ) 

+ + c 

so  that  'corrections'  satisfy  the  linear  equation 

V^6V  + n g"‘\6V  - 6V^)  = -V^'2  - 6"  (3.2) 

and  corresponding  (linearized)  boundary  conditions  on  fixed  boundaries. 

Equation  (3.2)  is  of  the  general  form 


Ll 


Lu  + Qu  = F (3.3) 


where  the  solution  u is  a function  of  x and  the  polar  angles  of  a point  on  the  unit  sphere, 
u = u(x,  P),  the  operator  L consists  of  two  parts 


_ 9 2 a 

with  ^ X ^ ^ ~ f = F(x,P)  are  known  functions.  Since  the  outside  bound- 

ax 

ary  conditions  take  the  form  (2.13),  a numerical  solution  of  the  problem  (3.  3)  must  be  sought  as  a 

truncated  expansion  (2.12),  i.e. 

N V N 

u (x,F)  = I C (X)  4>  (P)  , (3.-1) 

k=0 

where  N is  a suitable  cutoff  and  <tj,(P)  stands  (apart  from  a factor)  for  a (general)  spherical  harmon- 
ic Y^^(p,?>)  with  1-tm  being  an  even  number;  it  is  assumed  that 

Lp  $j^(P)  = $^(P)  (3.5) 

where  specializes  to  -f(f-t-l)  in  the  case  of  eg.  (3.2).  To  formulate  ' perturbed' Galerkin  methods 
for  obtaining  coefficients  C*^(x)  in  (3.4),  it  will  be  assumed  that  the  sequence  of  functions 
l^( P)  build  an  orthonormal  system  with  the  'discrete'  inner  product 

M(N) 

- |3-»> 

where  M(N)  points  P^  of  the  unit  sphere  and  weights  are  suitably  chosen. 

Existenceof  an  inner  product  of  fhe  form  (3.6)  for  a sequence  of  spherical  harmonics. 
{Y*^^(p,ip),  f-hm  an  even  number  , can  be  established  in  the  following  way.  It  is  known 

(see,  for  example.  Fox  and  Parker,  1968)  that  the  cosines  {cos  m build  an  orthogonal  system, 

because 


^ w cos  m <p,  cos  m'<p,  = f 
1=0  ’ ^ ^1 


7T 

with  <p.  = -j  and  w^  = 1 for  1 < 1 < J-1,  w^  = t = w^ 
see  Dahlquist  and  Bjorck  (19V  4).  Therefore, 


m # m' 


0 < m=m'  < J 


m = m'  = 0,J 


For  another  possible  choice  of  (w  ,<n  I-* 

1’  11=0 


I ij  0 < m < J 
I J m = 0,J 
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(J 


2f+l 

2ir  (/+m)'. 


J- 


4tt 


m ^ 0 


■ m = 0 


Since  the  product  p"^,(p)  is  an  even  polynomial  in  p of  degree  f+f ' which  is  even  and  at 


the  most  equal  to  2J,  there  is  an  exact  Gaussian  quadrature  such  that 

M, 

2 (f+m): 

2f+l  (f-m): 


I 


1=1  -1 


(Mj  is  the  minimal  number  of  points  for  which  the  above  relationship  is  true).  This  shows  that 


M 


J 


.5  )?0 


m 


where 


2tt 

m J 

J 

TT 

m = J 

J 

A = / 
m \ 


indicating  that  the  symbol  (P)  actually  stands  for  \^A  ^ ( p , <f ) The  number  of  discrete 

. , , K m f m 

(+) 

points  (p.,v>.)._|  always  exceeds  the  nujiiber  of  functions  in  the  sequence  {Y  (p, <,■’),  m + f 
J i 

an  even  integer}^  ^ for  a given  cutoff  J;  see  Table  1. 

Coefficients  of  the  'perturbed'  Galerkin  approximation  of  the  form  (3.  4)  to  the  solution  of  (3. 3) 
are  obtained  by  substituting  (3.  4)  into(3,  3)  and  taking  'discrete'  inner  products  of  the  resulting 


equation  and  ^’|^(P)  for  k = 0,1,.  . . ,N.  This  gives  the  following  system  of  N+1  equations  for  the 
N N 

coefficients  (x) 


(L  V Q*;  (X)C^X)=  F%, 


X 2 
X 


ks' 

s=0 


where 


M(N) 

. X w,  r(x,p,)*^(pp 


and  similarly. 


f =1 


2 % Q(x,P^)  $3(P,  ) <S>^  (P^) 


t =1 


(3.7) 


(3.8) 


(3.9) 


For  the  exact  (i.  e.  'unperturbed' ) Galerkin  method  the  matrix  and  the  vector  F^  are  defined  as 
corresponding  integrals  over  the  sphere  and  may  differ  In  value  from  the  expressions  given  by  (3.9) 
and  (3.8).  The  difference  can  be  displayed  by  expanding  the  functions  F and  Q 
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then 


Table  1.  is  the  total  number  of  coefficients  in  the  expansion  (2.  12)  with 
f < J.  M(J)  is  the  minimal  number  of  quadrature  points  in  the 
'discrete'  inner  product;  MU)  = 
oo 

F(x,P)  = V f (X)  $ (P)  , 
j=0  ^ ^ 

00 

Q(x,P)  = Yj  4>.(P)  , 

j=0  ^ ^ 


00  M(N) 

1=0  ^ f=l  ' J ' 


(P,) 


and 


M(N) 


oL'"'  ■ "<  ’k'h’nih 


) . 


(3.10) 


(3.11) 


,N 


The  value  of  vector  F , as  computed  by  (3.8)  or  (3.10),  can  be  made  arbitrarily  close  to  its  exact 


value  by  choosing  the  cutoff  N in  such  a way  that  f^(x)  for  j > N is  sufficiently  close  to  zero. 

N 

However,  even  if  for  this  choice  of  N , q^(x)  = 0 when  j > N,  tlie  matrix  Q as  computed  by  (3.9) 


or  (3,11)  will  contain  some  'aliasing'  terms. 

This  can  be  shown  by  assuming  that  the  'coordinate'  functions  3>.  (P)  couple  in  the  following 


*,(P) « (P) 
J m 


(3.12) 


j+m 

E 

s=  |j-m| 


C,  * (P) 
jms  s 


where 


Existence  of  such  a coupling  rule  for  the  spherical  harmonics  +m  even}  can  be  easily 

demonstrated  (see,  M.  E.  Rose,  1957),  Then,  denoting  by 


I^N) 


,^s  = E^ 


(3.13) 


N 

the  matrix  Q is  given  by 


N 


= E 


j=0 


)ks 


where  using  ( 3. 12) 


and 


^ 2N  M(N) 

iks  ■ Jio 


2N 


(3.14) 


In  other  words,  the  expressions  (3.13)  coincide  with  when  0 < j < N and  0 < k+s  < N,  but 

when  k+s  > N there  are  extra  terms.  Practice  shows  (see,  also,  Orszag,  1974)  that  these  aliasing 
terms  are  not  a serious  error.  It  is  always  possible  to  choose  M(N)  so  large  that  (3.6)  will  be  true 
for  0<k,s<2N  which  would  guarantee  “ ^jl^s  0<j,k,s<N.  In  that  case  equations  ( 3. 7 ) 

would  not  be  different  from  the  exact  Galerkin  approximation  to  (3.  3);  methods  of  this  type,  known  as 
spectral  methods,  have  been  successfully  applied  to  numerical  weather  prediction  (Bourke,  197  2)  and 
other  hydrodynamical  fluid  flow  problems  (Orszag,  197  4). 

Another  useful  form  of  'perturbed'  Galerkin  equations  is  obtained  by  multiplying  (3.7)  with 

a>  (P  ) and  summing  over  k ; the  result  is 
k m 


« X 2 

k=0  X 


(3.15) 


using  the  fact  that 
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The  second  term  in  (3.15)  can  be  written  as 


M(N)  N 


giving,  finally,  the  following  form  to  (3.15) 


X s=  1 


(3.16) 


B 


N 


+ Q(x,P^)u^x,P^)  = F(x,P^)  , 


N 


= )’  W (P  ) $,  (P  ) . 

ms  s k k s k m 


(3.17) 


In  this  formulation  the  method  is  known  as  the  orthogonal  collocation  method  (see,  Finlayson,  1972; 
Miketinac,  1976)  and  its  main  advantage  over  (3.7)  is  that  the  matrix  b'^  does  not  depend  on  x: 
the  disadvantage  is  that  the  number  of  unknowns,  i.e.  {u^(x,P  ) is  increased  to  M(N) 

Equations  (3.16)  can  be  derived  in  a different  way  starting  with  the  standard  collocation  method. 


The  question  of  convergence  of  methods  of  the  type  (3.7)  is  not  considered  here,  because 
there  is  already  an  extensive  amount  of  information  about  it  through  the  work  of  G.  M.  Vainikko  (1972, 
1967).  The  question  is  also  considered  by  Orszag  (1974). 

4.  Numerical  Procedure 

The  solution  of  structure  equations  is,  therefore,  sought  in  the  form 

V^x,p,^)  = j]  I 

/=0  m=0 
(/  +m  even) 

and  as  a limit  of  an  iterative  process  each  step  of  which  is  a two-p>oint  boundary  value  problem, 
consisting  of  a system  of  linear  ordinary  differential  equations.  For  the  'perturbed'  Galerkin  method 
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the  unknowns  are 


6vIq(0)  =n/2T  6 V ,{6v{  (X),  0 < X < X , /+m  even}J  ; 

- n ^f=0,m=0 


the  equations  are  for  0 < x < x 

2 ^ 
r d 2 d f(f+l)  , J 

s;  - -^1  * 


where 


(f  '+m'  even) 

Mj  j 

fJ  (X)  = - (ii+  2 i(l±l),vJ 

fm'  I 2^  X dx  2 J 

dx  X 


- L L e (x,ix 
i=l  j=0  ^ J ^ ^ J 


' Co-i-n'  • 


The  boundary  conditions  are,  in  their  linearized  form, 


d6vJo(0)_  dV„J„(0) 


- ■ 3’™'°'.  < ' » 


. 1±1  ...I  .,,  , “/m'V  „1,, 

H H 

An  equation  determining  6V^  is  obtained  by  integrating  (2. 1')  over  the  sphere  and  then  letting  x-  0. 
Since 

/d^/  (^(1-^2^  A+  — 1—  _!_jv(x,f,v>)  = 0 

0 I-|i  9^ 


''•’V'(x,„„  . ^ vj„(x,  . 


the  resulting  equation  is 


= - ’Jz]  , 


i 


-1^ 


which  in  the  linearized  form  is 

3 = -^^^T  - i y — . 

dx^  dx*^ 

For  the  method  of  orthogonal  collocation  the  unknowns  are 

M , J 

the  equations  are  for  0 < x < x 


dx  X p=l  g=l 


(4.7') 


+ n[e^(x,p.,?>j)]'’"  [6V(x,(i,,?>^)  - 6V^]  = R(x,p^,V>j)  , 

where 

= f,  y,  (-)f(l+l)A  w w ,v>  ) 

f=0  m=^0  ' m p q /m'^^p’^^q'  fm''^i'^j'’ 

(m+f  even) 

R(x,p.,?>J  H .(  + e"  ) 

+ V=Pj,»>=‘Pj 

The  boundary  conditions  are,  in  their  unlinearized  form, 

J . avVo,Pj,v.) 

L Tj  03^ = 0 . 

i=i  j=o  ^ 

1 i T T 

^ p^l  q-0  ‘ V^’^q  ’ = ° ' 

where 

S'j.pq  = ,i  i„  V'*"  "p  ’','1^-1  "i  I V',’- 


A (f+1)  w w y|'’'\p,  ,p,  ) Y^/^p  ,<p  ). 
m p q Im  i j fm  p’  q 


(4.10) 


(1.11) 


(4.12) 


(4.13) 


(4.14) 


To  determine  one  more  equation  must  be  supplied;  it  is  given  as  another  form  of  equation  (4.7) 


z’  i 

1 _i  i _n  * J 


dV(0,  p^,,,^) 


= -2J  . 


(4.15) 


1=1  j=0 


These  two  boundary  value  problems  are  singular,  but  there  is  a variety  of  methods  forobtaining 
their  numerical  so’ution  (see,  de  Hoog  and  Weiss,  1975;  also,  Barter,  1965).  Perhaps,  the  simplest 
method  is  to  use  central  finite  differences  on  an  equidistant  grid  x^  = h Ax,  h = 0,1,2, ...,  H:  this 
approach  is  adopted  in  this  investigation.  To  Improve  the  accuracy  differential  operators  involving 


the  corrections  (4.1)  on  the  left-hand  side  of  equations  (4.2),  (4.7'),  (4.5)  and  (4.6)  are  written  out 
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usirt^  the  simplest  central  finite  difference  formulas,  while  the  same  differential  operators  involving 
the  previous  iterate 

=1,2,...,H,  / + m even}/^’^^^^ 

on  the  right-hand  side  of  these  equations  are  represented  by  more  accurate  approximations  in  terms 
of  central  differences.  This  amounts  to  a fusion  of  Newton  corrections  (4.1)  and  deferred  corrections 
into  a single  iterative  procedure.  In  this  way  one  obtains  the  following  system  of  simultaneous  linear 
algebraic  equations 

and  for  h = 1,  2, . . . , H-1 


* l‘  ♦ ITl  t t J-hl- 

/ ' sO  m' =0  ' 


J 

z z 

!'=0  m'=C 
(f'+m'  even) 


(4.17) 


the  final  equation  is 


H 


H 


where 


__  ^ ±M  .r 

fm  dx  X tm  H 

H 


(4.18) 


(4.19) 


The  right-hand  sides  are  derived  in  the  appendix.  A similar  process  can  be  used  to  obtain  a numerical 
solution  of  equations  (4. 9),  (4.12),  (4. 13)  and  (4. 15). 

The  coefficient  matrix  for  the  algebraic  system  is  block-tridiagonal  except  for  the  full  first 
column  when  natural  ordering  of  unknowns  is  used  (see  Figure  1).  Such  systems  can  be  inverted 
directly  by  either  treating  the  coefficient  matrix  as  a band  matrix  or  by  exploiting  the  block-tridiag- 
onal structure.  However,  the  standard  methods  such  as  those  described  by  Martin  and  Wilkinson 


(1967)  or  Isaacson  and  Keller  (1966)  must  be  modified  to  take  into  account  the  full  first  column.  These 
modifications  are  particularly  simple,  if  the  ordering  of  unknowns  is  reversed.  Both  algorithms  have 


been  used  In  actual  computations  on  the  computer  taking  about  the  same  amount  of  computing  time  to 
Invert  a given  coefficient  matrix;  the  results  have  been  practically  identical  with  the  Issacson -Keller 
algorithm  being  more  efficient  because  the  amount  of  data  necessary  to  transfer  to  and  from  the 
massive  storage  is  halved.  In  spite  of  this  the  Martin-Wilkinson  algorithm  is  to  be  preferred  be- 
cause it  has  not  been  possible  to  show  that  the  condition  guaranteeing  numerical  stability  of  the 
Isaacson-Keller  algorithm  (see,  Varah,  197  2)  is  satisfied;  in  fact,  a simple  application  of  the 
Gerschgorin  theorem  shows  that  the  condition  is  not  satisfied,  but  it  is  very  close  to  it. 


L-  “h  *hJ 

Figure  1.  The  coefficient  matrix  for  equations  (4.16),  (4.17) 
and  (4.18).  All  blocks  are  square  matrices  of  order 
Nj  (Aj  has  one  additional  row  and  column)  except 
and  Cj.  The  order  of  the  coefficient  matrix  is 
H • Nj+ 1 ; c^ , c^, . . . , Cj^  j form  the  column  vector  in 
(4.17).  Blocks  Bj  (1  > 2)  and  (i  >1)  are  diagonal 
matrices. 

When  considering  the  computer  implementation  of  the  two  methods  described  above,  it  is  nec- 
essary to  take  into  account  the  cost  of  the  following  major  operations;  (i)  setting  up  the  coefficient 
matrix,  (ii)  solving  the  linear  system,  and  (iii)  storing  and  retrieving  data  from  the  massive  storage. 


There  are  H • 2(Nj)  Mj(J+l)  multiplications  involved  in  setting  up  the  coefficient  matrix  for  the 

3 

'perturbed'  Galerkin  method  and  H(Nj)  multiplications  needed  to  solve  the  corresponding  linear 
system  using  the  Martin-Wilkinson  algorithm;  the  total  is 
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(4.20) 


H(Nj)^(2Mj(J+l)  + Nj)  . 

For  the  method  of  orthogonal  collocation  the  number  of  multiplications  is  of  the  order  of  magnitudo 

H[Mj(J+l)]^  (4.21) 

since  setting  up  the  coefficient  matrix  is  not  a major  contribution  in  this  case.  Although  the  num;  i ; 
(4. 21)  is  greater  than  (4.  20)  for  J = 7 and  J = 10,  which  were  actually  used  on  the  computer,  the 
method  of  orthogonal  collocation  should  not  be  used  for  stellar  structure  computations  in  three  dimen- 
sions, because  the  cost  of  the  third  major  operation  mentioned  above  becomes  prohibitively  expen- 
sive very  rapidly  with  increasing  J . Also,  as  J increases  it  becomes  necessary  to  use  double 
precisioning  when  solving  the  linear  system  much  sooner  for  the  method  of  orthogonal  collocation 
bringing  a further  setbaclt  to  the  method. 

Even  for  two-dimensional  calculations  (Stoeckly,  1965),  where  the  number  of  coefficients 
equals  the  number  of  quadrature  points,  the  'perturbed'  Galerlcin  method  is  an  almost  equally  effi- 
cient alternative  to  the  method  of  orthogonal  collocation.  This  is  because  the  cost  of  the  major  oper- 
ations (i)  and  (ii)  is  about  equal  (the  factor  is  about  1.  5 if  symmetry  of  the  matrix  is  taken  into 

account),  while  the  cost  of  the  third  major  operation  does  not  have  to  be  considered  (even  thouglr  it 
is  the  same). 

5.  Results  of  Numerical  Experiments 

A computer  implementation  of  the  method  of  orthogonal  collocation  and  the  'perturbed'  Galerkln 
method  has  been  developed  for  the  problem  of  uniform  rotation  with  the  disturbing  function  given  by 
(2.7).  This  is  very  convenient  for  two  reasons.  The  rotating  configurations  should  be  axially  sym- 
metric, because  of  the  form  of  (at  least  for  small  £7,  see  James  , 1964),  so  that  the  problem  is 
actually  two-dimensional  and  as  such  its  numerical  solution  is  well-known.  So  the  present  methods 
will  be  thoroughly  tested  in  providing  a close  agreement  with  the  known  solution.  And,  finally,  the 
treatment  of  the  binary  problem  using  the  same  program  is  equivalent  to  a simple  replacement  of  the 
formula  (2.7)  by  (2.9)  or  (2,10)  in  addition  to  some  other  (minor)  complications.  However,  for  reasons 
already  mentioned  only  t)ie  'perturbed'  Galerkin  method  has  been  used  for  actual  computations  in  the 


binary  problem.  In  fact,  the  method  of  orthogonal  collocation  has  been  programmed  for,  mostly, 
sentimental  reasons  generated  by  its  superb  performance  in  two  dimensions  (Mlketinac,  1976),  due  in 


part  to  the  fact  that  there  the  number  of  coefficients  in  the  expansion  (3.4)  equals  the  number  of 
quadrature  (or  collocation)  points  in  (3.  6).  For  both  methods  the  basic  computations  can  be  organized 
as  in  the  flow  chart  in  Figure  2,  The  possibility  of  using  Ax  as  the  parameter  p in  the  flow  chart  | 

can  be  exploited  in  two  ways:  (i)  the  mesh  could  be  refined  (with  all  other  parameters  kept  fixed)  1 

j 

giving  an  idea  about  the  accuracy  of  an  obtained  solution,  and  (ii)  the  mesh  could  be  coarsened  to  1 

ensure  thot  the  maximal  radius  of  a particular  solution  is  smaller  than  H A x.  jt 

All  numerical  experiments  reported  here  have  been  performed  for  the  polytropic  index  n = 3.  0, 
with  H equal  to  60,  and  J has  been  chosen  to  be  equal  to  7 except  for  a few  experiments  when  } | 

has  been  equal  to  10.  Table  2 contains  some  pertinent  information  about  the  storage  requirements  j 

and  computing  times  involved  in  an  iteration  cycle  for  the  'perturbed'  Galerkin  method.  In  the  case  ' 

of  the  method  of  orthogonal  collocation  it  has  been  necessary  to  use  double  precisioning  even  for  J =7;  | 

the  computing  time  per  iteration  cycle  has  been  about  I minute.  From  the  Table  2 it  is  clear  that  fre- 
quent use  of  the  massive  storage  must  be  made.  This  is  achieved  by  modifying  the  Ma.nin -Wilkinson  | 

algorithm  (for  which  the  coefficient  matrix  is  a band  matrix  of  band  width  2N  +1  and  having  60-  N +1  |j 

J J ij 

rows)  in  such  a way  that:  (i)  the  coefficient  matrix  is  set  up  in  consecutive  blocks  containing  Nj  |i 

rows  each  (except  the  first  block  which  contains  N^+l  rows),  (ii)  the  fast  core  contains  only  two 
such  blocks  at  any  one  time  and  (iii)  results  of  the  elimination  are  transferred  to  and  from  the  massive 
storage  in  blocks,  too.  For  all  values  of  the  parameter  p in  Figure  2 that  have  been  tried,  the  num- 
ber of  iterations  necessary  for  convergence  (convergence  criteria  are  that  max  IR  (x  )1  <0  000  001 

f,m,h  fm  h ' - 

and  max  l6V,,^,(x.)l  < 0.000  6)  seldom  exceeded  three. 

/,m,h  fM  h ' - 

Computation  of  uniformly  rotating  configurations  has  served  (as  indicated  at  the  beginning)  as 
a test  problem  in  two  ways.  The  accuracy  of  a method  can  be  judged  by  (i)  looking  for  dependence  ] 

J 

of  the  results  on  the  angle  <p  , and  (it)  agreement  with  the  corresponding  two-dimensional  results.  ' 

For  the  'perturbed'  Galerkin  method  no  dependence  on  <p  is  found  to  seven  decimal  digits  and  the 
results  agreed  with  the  two-dimensional  results  (obtained  by  the  method  of  orthogonal  collocation 
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I 

"s 

"r 

T 

7 

20 

50  360 

32 

8 

40 

10 

36 

159  768 


158 



29 

187 

Table  2.  The  new  symbols  l^j3>Tg,Tj^,  and  T have  the  following  meanings: 

Nj^  is  the  number  of  non-zero  elements  for  the  Martin-Wilkinson 

algorithm,  Tg  is  the  time  taken  for  setting  up  the  coefficient 

matrix  and  solving  for  the  corrections,  T is  the  time  taken  for 

R 

calculating  the  right-hand  sides,  and  T is  the  total  time  per 
iteration  cycle;  all  times  are  in  seconds  (the  machine  is  UMVAC  1110). 

formulated  for  two-dimensional  computations,  Miketinac  and  Barton,  19~2)  to  four  decimal  digits. 

These  results  have  been  slightly  bettered  by  the  method  of  orthogonal  collocation  out  with  the  Gav.ss 

elimination  done  in  double  precision.  Both  computations  have  been  performed  for  J = T made  possible 

by  the  fact  that  the  magnitude  of  the  coefficients  V,  (x.  ) becomes  very  small  for  f >S.  Marti:. 

< m h — 

(1970),  also,  tested  his  perturbation  expansion  on  the  example  of  uniformly  rotating  polytropos.  His 
results  for  n = 3.0  agreed  very  well  with  the  results  obtained  by  Mikctinac  and  Barton  (1972)  usr.c 
the  method  of  orthogonal  collocation:  the  results  are  shown  in  Table  3. 


-3 

- ^ 

n 

■1.  070X  10 

4. 076  X 10 

c 

X 

10.05 

10.065 

e 

Table  3.  Equatorial  radius,  , for  the  polytrope  of  index 
n = 3.0  rotating  uniformly  with  the  critical  angular 
velocity,  The  second  column  contains  results 

obtained  by  Mikctinac  and  Barton,  while  the  third 
column  contains  Martin's  results. 
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Computation  of  the  structure  of  the  primary  in  a close  binary  system  has  been  perforir.ed,  so 
far,  for  only  one  value  of  the  ratio  of  masses,  q = 0.  5 (with  n = 3.0).  The  first-order  theory  has 
been  used,  J being  equal  to  7 and  the  critical  separation  is  found  by  testing  (2.U).  The  re- 
sults are  collected  in  the  second  column  of  Table  4,  which  contains  (in  the  third  column)  the  corrc - 
spending  results  obtained  by  P.  G.  Martin  (1970)  using  the  second-order  theory  and  the  perturbatio- 
method.  The  fact  that  results  do  not  agree  much  more  closely  can  probably  be  accounted  for  by  the 
difference  between  the  first-order  theory  and  the  second-order  theory.  Tor  the  critical  configuratio  . 

the  absolute  magnitude  of  coefficients  V.  (x,  ) becomes  smaller  than  0.000  002  for  f >6  shov.ur.c 

fm  h — ■ 

that  the  cutoff  J = 7 is  reasonable.  An  intermediate  configuration,  belonging  to  .\  * = 0.03,  has 
been  computed,  also,  with  J = 10  (and  still  using  the  first-order  theory),  the  results  agreed  to  five 
decimal  digits.  Another  intermediate  configuration,  belonging  to  X * = 0.06,  has  been  obtained  in 
two  different  ways:  (i)  as  an  end  of  a sequence  of  configurations  with  X ^ making  'small  jumps' 

0.  00  — 0.  03  -♦  0.  04  -♦  0.  05-*  0.  06,  and  (ii)  directly  from  X ^ = 0.  00  to  X ' = 0.  Ou;  the  results 
agreed  again  to  five  decimal  digits. 

The  immediate  future  work  on  the  binary  problem  should  have  the  following  two  objectives. 

First,  critical  configurations  for  a range  of  values  of  the  mass  ratio,  q (including  q = 0.  h)  slioulv 

be  obtained  using  both  the  first-order  and  the  second-order  theory  with  J = 10  and  the  results  shoul  ' 

then  be  compared  v.ith  those  obtained  by  Martin  (19701:  a close  agreement  of  the  results  is  expected. 

In  these  computations  the  polytropic  index  would  have  to  be  equal  to  n = 3.0,  since  Martin's  method 

can  not  be  extended  to  the  other  interesting  polytrope  for  v,’hich  n = 1.5.  The  second  objective  sho-ld 

be  computation  of  critical  configurations  for  this  other  polytropc  using  again  both  the  first-order  and 

the  second-order  theory  with  J = 10.  The  results  could  then  be  compared  with  the  first-order  results 

of  Naylor  and  Anand  (1970),  but  large  discrepancy  of  the  results  would  not  be  surprising  on  account 

of  remarks  made  in  Section  2.  It  is  expected  that  H = 60  will  have  to  be  replaced  by,  probably, 

H = 100  in  order  to  maintain  A x < 0.  2 and  yet  accommodate  the  configurations  for  the  n = 1.  5 

polytrope,  which  is  known  to  be  less  centrally  condensed  than  the  n = 3.0  polytrope.  This  change 

In  II  will  make  the  computations  more  expensive  and  an  increased  efficiency  will  have  to  be 

achieved.  For  example,  the  cost  of  com.puting  critical  configurations  using  the  second -order  theory 

-21- 


x-^ 

C 

0.  064 

0.  062 

X 

e 

0.  583 

0.  57  1 

0.  459 

0.  439 

X 

n 

0.  497 

0.  468 

X 

P 

0.  428 

0.414 

"p 

1.092 

1.  09  4 

^t 

0.756 

0.  780 

^n 

0.  533 

0.  601 

V 

s 

2.  291 

2.  961 

Table  4.  A comparison  of  results  for  the  primary  of  a close  binary 

system.  The  symbols  x ,x  ,x  , and  x are  the  radii 

e t n p 

in  the  following  (p,  i?)-directions : (0,0),  (0,tt/2).  (0,’t), 

and  {l,<p).  The  gradient  of  the  (total)  potential  (i.  e.  V+V^^) 

at  the  surface  is  denoted  by  g;  the  quantities  g and  g 

t 'n 

are  expressed  in  the  units  of  g and  g is  itself  normal- 

P P 

ized  by  the  gravity,  g^,  on  the  undistorted  Emdcn  sphere. 
The  (total)  potential  at  the  surface  is  the  potential 
normalized  by  Gm/R. 


with  J = 10  can  be  lowered  considerably  by  first  locating  the  critical  configuration  in  the  first-order 
theory'  with  J = 7 and  then  refining  the  results  using  the  second-order  theory  and  J = 10. 

The  future  work  on  the  'perturbed'  Galerkin  m.ethod  should  include  the  proof  of  convergence. 
Cryer  (1976)  reports  that  no  proof  of  convergence  of  a trial  free-boundary  method  has  been  given. 

From  the  applications  point  of  view  the  most  important  future  problems  would  be  the  bifurcation  prob- 
lem for  the  disturbing  potential  (2.8)  and  the  possibility  of  computing  the  critical  configurations  of 
more  realistic  models  of  stars. 
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Appendix 


From  equations  (-1.  ?>)  and  (-1.7')  it  follows  that  the  right-hand  side  in  (-4. 16)  is  given  by 

■■oV”'  h 

dx 

The  right-hand  side  for  (4.17)  is  defined  in  (4.4),  and  the  final  right-hand  side  in  (4.18)  is  a combina- 
tion of (4.4) and  (4.19). 

Differential  operators  in  these  right-hand  sides  are  approximated  using  the  following  formulas 


12(Ax) 


r<-V2-'‘^yi+i  - '°yi  + ^S-i  -^-2*- 


Then  one  can  show  that  the  right-hand  sides  are  given  b>' 

’’L""  - . ,0  V„'„(0)  t , 

> TSlii*  K'  ■ '<■"  ‘ 

+ 30  v/  (X  ) - 16(1  -^)v/  (x^  ,)+  (1  - ^)V/  (x^  \]  + 

rmh  hfmh-1  hfmh-2 

] 

^ '"i^j  °>h'^’^i’ 

for  h = 1, 2, . . . , H-1 , and 

* ""*  5'“  °,m  - II  * 

H 

These  expressions  can  be  used  only  if  a prescription  is  given  for  computing  V^^fx  ^),  V^^(x  j)  and 
^ quantities  are  obtained  from  the  outside  boundary  condition  by 


setting  In 


^Im^’^M+l^  ■ ^fm^’^M-l'  f+1  .-J  . 

2^^; ^1^  ''lm<’'M>  = ° 


M = n and  M = H+1.  An  expression  for  V (x  ) is  obtained  by  putting  h = -1  in  the  equation 

f m -2 
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J. 


' W»t,>  - “ -h-i  ' 


where 


7 J 


i=l  j=0 


i j 00''^i’'j'  +' 


giving 


0(x_j,fi,v>)  = V(x_j,^,v>)  - + 1 + Vj^(\,x_j,^,v>)  , 


2 


It  follows  from  (4.1)  that  when  1=0  one  can  set 

v/„(x  ,)  = (-)^  ■ 

tm  -l  f m 1 

The  validity  of  this  formula  for  1 0 follows  from  the  known  solution  (Weinberger,  196  5)  of 

d^V  d V 

Im  2 fm  f (f +1)  - _ 

.2  xdx  2 fm" 

dx  X 

near  x = 0. 
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